library(tidyverse)
library(marginaleffects)
library(ggsci)
library(Cairo)
library(GPArotation)
library(estimatr)
library(interactions)
library(mediation)
library(tibble)
library(Hmisc) 

#csvデータを読み込む
data <- read_csv("dataset_sekai.csv")

data <- data %>% 
  mutate(sex = case_when(Q3.1 == 1 ~ "女性",
                         Q3.1 == 2 ~ "男性",
                         TRUE ~ NA)) %>% 
  mutate(pid = 
           case_when(Q9.1 == 1 ~ "自民",
                     Q9.1 == 2 ~ "立憲",
                     Q9.1 == 3 ~ "維新",
                     Q9.1 == 4 ~ "公明",
                     Q9.1 == 5 ~ "国民",
                     Q9.1 == 6 ~ "共産",
                     Q9.1 == 7 ~ "社民",
                     Q9.1 == 8 ~ "れ新",
                     Q9.1 == 9 ~ "参政",
                     Q9.1 == 10 ~ "保守",
                     Q9.1 == 12 ~ "支持なし",
                     TRUE ~ NA),
         pid = factor(pid, levels = c("自民","立憲","国民","公明","維新","共産","社民",
                                      "れ新","参政","保守","支持なし")),
         財務省 = Q7,
         age = if_else(Q3.2 == 999, NA, Q3.2),
         edu = if_else(Q3.3 == 9, NA, Q3.3),
         income = if_else(Q3.4 == 99, NA, Q3.4),
         inte = if_else(Q4.1 == 9, NA,5 - Q4.1))

data <- data %>% 
  mutate(con1 = Q11.1_1,
         con2 = Q11.1_2,
         con3 = Q11.1_3,
         con4 = Q11.1_4,
         con5 = Q11.1_5) %>% 
  mutate(con = con1 + con2 + con3 + con4 + con5,
         conz = scale(con))


data <- data %>% 
  mutate(X = if_else(Q8.1_4 == 9, NA, 4 - Q8.1_4),
         youtube = if_else(Q8.1_5 == 9, NA, 4 - Q8.1_5),
         TikTok = if_else(Q8.1_6 == 9, NA, 4 - Q8.1_6)) 

########################################################################################
#Figure1

pivot_a <- data %>%
  filter(!is.na(財務省), !is.na(we)) %>%
  pivot_longer(cols = c(財務省),
               names_to = "group", values_to = "value") %>%
  group_by(pid, group) %>%
  summarise(
    mean = weighted.mean(value, w = we),
    var = Hmisc::wtd.var(value, weights = we),
    n = n(),
    sd = sqrt(var),
    se = sqrt(var / n),
    ci = se * 1.96,
    .groups = "drop"
  )


mean_weighted1 <- data %>%
  filter(pid != "国民" | pid != "れ新" | pid != "参政") %>% 
  drop_na(財務省,we) %>%
  summarise(mean = weighted.mean(財務省, w = we)) %>%
  pull(mean)
mean_weighted1

mean_weighted <- data %>%
  drop_na(財務省,we) %>%
  summarise(mean = weighted.mean(財務省, w = we)) %>%
  pull(mean)
mean_weighted

pivot_a　<- pivot_a %>% 
  mutate(check = if_else(pid == "国民" |pid == "れ新" |pid == "参政", "sig", "ns"))

Fig1 <- pivot_a %>% 
  drop_na(pid) %>% 
  ggplot(aes(x = pid, y = mean, color = check)) +
  geom_pointrange(aes(ymin = mean - ci, ymax = mean + ci)) +
  labs(x = "支持政党", y = "財務省好感度の平均値 w/ 95%信頼区間")+
  theme_bw(base_size = 12, base_family = "HiraginoSans-W3") +
  geom_text(aes(x = pid, y = 1, label=sprintf("%.2f", mean)),
            size = 4, color = "black") + 
  guides(color = "none") + 
  geom_hline(yintercept = mean_weighted, linetype = "dashed", color = "black") + 
  theme(axis.title.y = element_text(size = 10),  
        axis.title.x = element_text(size = 11)) + 
  scale_color_manual(values = c("sig" = "black",
                                "ns" = "grey70"))
  

Fig1
ggsave("Figure1.png", Fig1, dpi = 300, width = 9, height = 3)

################################################################################################
#Figure2

data2 <- data %>%
  dplyr::select(財務省, conz, pid, sex, age, edu, income, inte, we) %>%
  na.omit() %>%
  mutate(pid = factor(pid))  # pidはfactorのままでOK

# 対象の政党
parties <- c("立憲","国民","維新","共産","社民","れ新","参政","保守")

# 結果を保存するリスト
results <- list()

for (party in parties) {
  data_tmp <- data2 %>%
    mutate(pid_bin = if_else(pid == party, 1, 0))  # ダミー変数
  
# モデル推定
model.m <- lm(conz ~ pid_bin + sex + age + edu + income + inte, data = data_tmp)
model.y <- lm(財務省 ~ pid_bin + conz + sex + age + edu + income + inte, data = data_tmp)
  
# 媒介分析
med.out <- mediate(model.m, model.y, treat = "pid_bin", mediator = "conz", sims = 1000)
  
# 結果抽出
summary_out <- summary(med.out)
results[[party]] <- tibble(
    pid = party,
    estimate = summary_out$d0,
    lower = summary_out$d0.ci[1],
    upper = summary_out$d0.ci[2],
    pval = summary_out$d0.p
  )
}



# データフレームにまとめる
df_prop <- bind_rows(results) %>%
  mutate(check = ifelse(pval < 0.05, "sig", "ns")) %>% 
  mutate(pid = factor(pid, levels = c("立憲","国民","維新","共産","社民","れ新","参政","保守")))

mean_weighted_con <- mean(df_prop$estimate, na.rm = TRUE)  # 中央線の参考

Figure2 <- df_prop %>%
  ggplot(aes(x = pid, y = estimate, color = check)) +
  geom_pointrange(aes(ymin = lower, ymax = upper)) +
  geom_text(aes(x = pid, y = estimate  , label = sprintf("%.2f", estimate)),
            size = 4, hjust = -.3, color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray70") +
  scale_color_manual(values = c("sig" = "black", "ns" = "grey70")) +
  labs(x = "支持政党", y = "平均因果媒介効果（ACME） w/ 95%信頼区間") +
  theme_bw(base_size = 14, base_family = "HiraginoSans-W3") +
  theme(
    axis.title.y = element_text(size = 11),
    axis.title.x = element_text(size = 12)
  ) + 
  guides(color = "none")


Figure2

ggsave("Figure2.png", Figure2, dpi = 300, width = 11, height = 4)

###############################################################################################################
#Figure3

data <- data %>% 
  drop_na(con1,con2,con3,con4,con5)

con_mean_weighted <- data %>%
  drop_na(conz,we) %>%
  summarise(mean = weighted.mean(conz, w = we)) %>%
  pull(mean)

con_mean_weighted

model <- lm_robust(
  財務省 ~ sex + age + edu + income + inte + 
    youtube * conz + TikTok * conz + X * conz,
  data = data,
  se_type = "HC2",
  weights = we)

slope_y <- slopes(
  model,
  variables = "youtube",
  newdata = datagrid(
    conz = seq(min(data$conz), max(data$conz), by = 0.05),
    model = model            # model を渡せば他の変数が自動で平均/最頻値に設定される
  ),
  conf_level = 0.95) %>% 
  mutate(group = "YouTube") %>% 
  dplyr::select(conz, estimate, conf.low, conf.high, group)

slope_tv <- slopes(
  model,
  variables = "TikTok",
  newdata = datagrid(
    conz = seq(min(data$conz), max(data$conz), by = 0.05),
    model = model            # model を渡せば他の変数が自動で平均/最頻値に設定される
  ),
  conf_level = 0.95) %>% 
  mutate(group = "TikTok") %>% 
  dplyr::select(conz, estimate, conf.low, conf.high, group)

slope_x <- slopes(
  model,
  variables = "X",
  newdata = datagrid(
    conz = seq(min(data$conz), max(data$conz), by = 0.05),
    model = model            # model を渡せば他の変数が自動で平均/最頻値に設定される
  ),
  conf_level = 0.95) %>% 
  mutate(group = "X（旧Twitter）") %>% 
  dplyr::select(conz, estimate, conf.low, conf.high, group)


slope <- rbind(slope_y, slope_tv, slope_x) %>% 
  mutate(group = factor(group, levels = c("YouTube","TikTok","X（旧Twitter）")))


Figure3 <- slope %>%
  ggplot(aes(x = conz, y = estimate)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              fill = "gray80", alpha = 0.7) +
  geom_line(color = "black", size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    x = "陰謀論志向",
    y = "各サービス接触が財務省好感度に与える影響"
  ) +
  theme_bw(base_size = 12) +
  facet_wrap(~ group) + 
  theme(axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 11))


Figure3
ggsave("Figure3.png", Figure3, dpi = 300, width = 8.5, height = 3)


##########################################################################################################################
# Figure4

model <- lm_robust(
  財務省 ~ sex + age + edu + income + inte + 
    youtube * conz + TikTok * conz + X * conz,
  data = data,
  se_type = "HC2",
  weights = we)

pred1a <- predictions(
  model,
  newdata = datagrid(
    conz = c(min(data$conz, na.rm = TRUE), max(data$conz, na.rm = TRUE)),
    youtube = c(min(data$youtube, na.rm = TRUE), max(data$youtube, na.rm = TRUE)),
    model = model )) %>% 
  mutate(group = "YouTube",
         media = youtube) %>% 
  dplyr::select(estimate,conz,media,conf.low,conf.high,group)

pred2a <- predictions(
  model,
  newdata = datagrid(
    conz = c(min(data$conz, na.rm = TRUE), max(data$conz, na.rm = TRUE)),
    X = c(min(data$X, na.rm = TRUE), max(data$X, na.rm = TRUE)),
    model = model )) %>% 
  mutate(group = "X（旧Twitter）",
         media = X) %>% 
  dplyr::select(estimate,conz,media,conf.low,conf.high,group)

pred3a <- predictions(
  model,
  newdata = datagrid(
    conz = c(min(data$conz, na.rm = TRUE), max(data$conz, na.rm = TRUE)),
    TikTok = c(min(data$TikTok, na.rm = TRUE), max(data$TikTok, na.rm = TRUE)),
    model = model )) %>% 
  mutate(group = "TikTok",
         media = TikTok) %>% 
  dplyr::select(estimate,conz,media,conf.low,conf.high,group)

pred <- rbind(pred1a,pred2a,pred3a) %>% 
  mutate(conz = factor(conz, labels = c("最小","最大")),
         media = factor(media, labels = c("全く利用しない","よく利用"))) %>% 
  mutate(group = factor(group, levels = c("YouTube","TikTok","X（旧Twitter）")))


Figure4 <- pred %>%
  ggplot(aes(x = conz, y = estimate, shape = media)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  width = 0.2,
                  position = position_dodge(width = 0.5)) +  
  labs(x = "陰謀論志向", y = "財務省好感度の予測値 w/ 95%信頼区間",
       shape = "メディア接触量") +
  theme_bw(base_size = 12)  + 
  theme(legend.position = "bottom",
        axis.title.y = element_text(size = 10),  
        axis.title.x = element_text(size = 11),
        legend.box.margin = margin(t = -10)) +
  facet_wrap(~ group) + 
  geom_text(pred　%>% filter(media == "全く利用しない"),
            mapping=aes(x = conz, y = estimate,
                        label = sprintf("%.2f",estimate)),
            size = 3, nudge_x = -.3, nudge_y = 0, segment.size = 0) +
  geom_text(pred　%>% filter(media == "よく利用"),
            mapping=aes(x = conz, y = estimate,
                        label = sprintf("%.2f",estimate)),
            size = 3, nudge_x = .3, nudge_y = 0, segment.size = 0) + 
  scale_y_continuous(breaks = seq(0, 80, by = 10),limits=c(0, 80))

Figure4
ggsave("Figure4.png", Figure4, dpi = 300, width = 10, height = 4)

###########################################################################

